home *** CD-ROM | disk | FTP | other *** search
/ Ian & Stuart's Australian Mac 1993 September / September 93.iso / Archives / Applications / Calculators / NumberCrunch 1.41 / Number Crunch II Demo / Writing External Codes (xCOD) / FORTRAN source code / PROPOGATE.f < prev    next >
Encoding:
Text File  |  1992-08-16  |  1.8 KB  |  68 lines  |  [TEXT/MPS ]
  1. C This can be compiled into an xCOD resource and inserted into NumberCrunch.
  2. C It needs an sCOD resource with the same ID, which will contain the Ascii characters:
  3. C
  4. C     xCOD xPROPOGATE(dT,zN:extended ;  X,Y: array ; program VXY(n:extended; xy:array) ), propogates a classical trajectory 
  5. C
  6. C
  7.  
  8.     SUBROUTINE PROPOGATE(dT, zN, X, Y, VXY)
  9.     implicit EXTENDED (a-h,o-z), INTEGER (i-n)
  10.     dimension X(*), Y(*)
  11.     external VXY
  12.  
  13. C PROPOGATE will send a particle forward along a 2 dimensional
  14. C classical trajectory. The gradient of the potential energy
  15. C (which is the negative of the force) is calculated by 
  16. C the external subroutine or user program VXY, 
  17. C which is called via
  18. C    Varray[1] = x
  19. C    Varray[2] = Y
  20. C    CALL VXY(2.0,Varray)
  21. C  Varray has (x,y) as input to VXY, and returns (dV/dX, dV/dY) as output.
  22. C 
  23. C Input to xPROPOAGATE:
  24. C            dT is the time step,
  25. C            zN is the size of the X and Y arrays,
  26. C            X(1), X(2), Y(1), and Y(2) are set to start the trajectory,
  27. C            VXY is an external subroutine or user program to calculate the gradient.
  28. C Output:
  29. C           X(3 to zN) and Yarr(3 to zN) have been calculated.
  30. C
  31. C Note that 
  32. C        - X and Y MUST point to zN component arrays (or the system crashes)
  33. C       - (zN-3) time steps forward are calculated.
  34. C
  35. C The algorithm used is just:
  36. C    X(n+1) = 2*X(n) - X(n-1) - dT^2 * Vx( X(n), Y(n) )
  37. C    Y(n+1) = 2*Y(n) - Y(n-1) - dT^2 * Vy( X(n), Y(n) )
  38. C
  39. C which you get from F=ma with m=1 , F=[-Vx,-Vy] , and 
  40. C    aX(n) = [X(n+1) - 2X(n) + X(n-1)]/dT^2, 
  41. C which is the simplest possible discretization.
  42. C
  43. C                        -Jim Mahoney,  7/11/89
  44.  
  45. C
  46. C  
  47.     extended Varray(2)
  48.     extended Two
  49.     
  50.     Two = 2.0
  51.     
  52.     Nmin1 = zN-1
  53.     dT2 = dT*dT
  54.     
  55.     do i=2,Nmin1
  56.     
  57.         Varray(1) = X(i)
  58.         Varray(2) = Y(i)
  59.         call VXY(Two, Varray)
  60.     
  61.         X(i+1) = 2.0*X(i) - X(i-1) - dT2 * Varray(1)
  62.         Y(i+1) = 2.0*Y(i) - Y(i-1) - dT2 * Varray(2)        
  63.     end do
  64.  
  65.     return
  66.     end
  67.  
  68. C That's it…